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ABSTRACT 

We investigate the atomic (HI) and molecular (H 2 ) Hydrogen content of normal 
galaxies by combining observational studies linking galaxy stellar and gas budgets to 
their host dark matter (DM) properties, with a physically grounded galaxy formation 
model. This enables us to analyse empirical relationships between the virial, stellar, 
and gaseous masses of galaxies and explore their physical origins. Utilising a semi- 
analytic model (SAM) to study the evolution of baryonic material within evolving DM 
halos, we study the effects of baryonic infall and various star formation and feedback 
mechanisms on the properties of formed galaxies using the most up-to-date physical 
recipes. We find that in order to significantly improve agreement with observations 
of low-mass galaxies we must suppress the infall of baryonic material and exploit a 
two-phase interstellar medium (ISM), where the ratio of HI to H 2 is determined by 
the galactic disk structure. Modifying the standard Schmidt-Kennicutt star formation 
law, which acts upon the total cold gas in galaxy discs and includes a critical density 
threshold, and employing a star formation law which correlates with the Hi gas mass 
results in a lower overall star formation rate. This in turn, allows us to simultaneously 
reproduce stellar, HI and Hi mass functions of normal galaxies. 
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1 INTRODUCTION 

Neutral atomic hydrogen is the most abundant element in 
the Universe and plays a fundamental role in galaxy for- 
mation, principally as the raw material from which stars 
form. Within galaxies, the Interstellar Medium (ISM) acts as 
a temporally evolving baryonic component; competing pro- 
cesses cause the accumulation (through external infall from 
the intergalactic medium and stellar evolution) and deple- 
tion (through star formation and various feedback mecha- 
nisms) of hydrogen. Thus, observational determinations and 
theoretical predictions of the hydrogen budget within galax- 
ies of various masses and morphologies is of central impor- 
tance to constraining the physics of galaxy formation (see 
Kauffmann, White & Guideroni, 1993, Benson et al. 2003, 
Yang, Mo & van den Bosch, 2003, Mo et al. 2005, Kaufmann 
et al. 2009) 
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Moreover, within the ISM Hydrogen comprises the ma- 
jority of the cold gas mass, and when non-ionized ex- 
ists within two-phases, atomic HI and molecular H2 ■ A 
large body of observational analysis has shown that within 
galaxies, HI generally follows a smooth, diffuse distribu- 
tion whereas H2 regions are typically dense, optically thick 
clouds which act as the birthplaces for newly formed stars 
(Drapatz & Zinnecker, 1984, Wong & Blitz, 2002, Krumholz 
& McKee, 2005, Blitz & Rosolowski, 2004, 2006, Wu et 
al., 2005). Due to the distinct differences in these phases, 
and the central importance of ISM physics to the evolution 
of galaxies, cosmological simulations have begun to include 
both phases (see Gnedin et al. 2009 & references therein), 
and observations have begun focusing on simultaneous mea- 
surements of both HI and H2 (see Obreschkow & Rawlins, 
2009). 

The distinction between these two phases has recently 
been shown to be of crucial importance to constrain the 
physics of galaxy formation, in particular resolved spec- 
troscopy using GALEX showing obscured star forming re- 



2 M. Cook et al. 



DM merger-trees 



Main Progenitor 
MAH 



Halo properties: 
2-phase evolution 



Fast accretion 
phase 



Spheroid-SMBH 
coevolution 



Slow accretion 
phase 



Disk growth 



z=0 Galaxy 
properties 



Figure 1. Schematic of the model framework, showing how we 
partition the evolutionary history into spheroid and disc growth 
epochs, which follow significantly different evolutionary paths for 
baryonic structure growth and result in the bulge-disk dichotomy 
observed at z = 0. 



gions in nearby galaxies (Kennicutt et al. 2003, 2007, 
Calzetti et al. 2007, Gil de Paz et al. 2007), and various 
observational surveys providing maps of gas in galaxies at 
high-resolutions (Walter et al. 2008, Heifer et al. 2003, Leroy 
et al. 2009), have revealed a deeper level of complexity on 
sub-galactic scales. These studies allowed theoretical mod- 
els for the ISM and star formation to be constrained and 
further developed. 

Furthermore, due to the constant replenishment and de- 
pletion of Hydrogen in either HI or H2 phases, and to their 
separate yet interlinked properties, at any epoch, measure- 
ments of the fraction of HI and H2 are highly constraining 
for the processes of molecular cloud formation, star forma- 
tion, baryonic infall and various feedbacks. Therefore, simul- 
taneously predicting the stellar and gas mass functions of 
normal galaxies is a major challenge for any physically mo- 
tivated galaxy formation model, requiring an accurate depic- 
tion of all of the aforementioned processes (see Mo et al. 2005 
for a detailed discussion) . These issues manifest most clearly 
within the largely successful ACDM paradigm within the 
lowest mass systems, where it still remains unclear whether 
strongly non-linear feedback mechanisms, lower star forma- 
tion efficiencies, or suppression of initial infall onto DM ha- 
los is the dominant driver for the suppression of luminous 
structure formation (Mo et al. 2005). It is more than likely 
that a combination of the above-mentioned effects will go a 



long way to alleviating current tensions between models and 
observations, since current semi analytical models (SAMs) 
incorporate several processes in order to generate a defi- 
ciency of stellar mass in DM halos; many of which operate 
most effectively at low masses (Benson et al. 2003, De Lucia 
et al., 2004). 

Observationally, Zwaan et al. (2005) used the catalogue 
of 4315 extragalactic HI 21-cm emission line detections from 
the HI Parkes All Sky Survey (HIPASS, Barnes et al. 2001) 
and obtained the most accurate measurement of the HI 
mass function of galaxies to date. The HIMF is fitted with 
a Schechter function with a faint-end slope of —1.37 ± 0.03. 
The sensitivity of this survey was so high that they were 
able to extend their analysis well down to HI masses of 
10 7 ' 2 Mq, hence this is most complete analysis so far. Using 
these statistical constraints, it has now become possible to 
make stringent comparisons between theoretical models and 
observations even in low mass galaxies. 

The physics of cold gas becomes increasingly relevant 
for constraining galaxy formation models at relatively low 
masses (dominated by late- type galaxies), where the pres- 
ence of gas becomes substantial and therefore may break 
the degeneracies between feedback, star formation, and in- 
fall processes. Moreover, within the ACDM scenario the HI 
and H2 mass budgets in galaxies are determined by an intri- 
cate offset between several competing processes, all of which 
have strong mass dependencies. Thus the present HI and H2 
fractions are strong functions of host DM halo mass and the 
evolutionary history of each individual galaxy. More specif- 
ically, the fraction of gas which may be captured by the 
host DM halo, and in turn removed by feedback, is expected 
to depend strongly on the binding energy of the gas itself, 
which is principally determined by the DM halo virial mass 
and density distribution. Thus, under this framework the 
observational properties of galaxy populations are strongly 
influenced by their collective host DM halos (White & Rees, 
1978, see Somerville et al. 2008 for a review). 

Motivated by the above-mentioned observational ad- 
vances and theoretical challenges, the primary aim of this 
work is to investigate the physical origins of the relation- 
ships between HI , stellar, and virial masses of galaxies. In 
order to do this we compare empirical galaxy relations de- 
rived from observational studies, with a physically motivated 
galaxy formation model. Observationally, we use a numeri- 
cal approach (described in Shankar et al. 2006) that relies 
on the assumption of the existence of a one-to-one mapping 
between galaxy properties and host DM halo mas^j. We 
interpret these results using a physically motivated SAM, 
(see Cook et al. 2009, hereafter C09, & references therein), 
which has been shown to reproduce many features of the 
local galaxy population. 

SAMs provide a powerful theoretical framework within 
which we can explore the range of physical processes (e.g. 
accretion mechanisms, star formation, SN feedback, black 
hole growth and feedback etc.) that drive the formation and 
evolution of galaxies and determine their observable proper- 
ties (see Somerville et al. 2008 and Baugh 2006 for extensive 



1 This approach is based also on the assumption on the com- 
pleteness of the sample over which the mass functions has been 
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reviews) whilst remaining computationally feasible, allowing 
for statistical samples of galaxies to be generated. 

As is shown in Fig.l, the backbone to our galaxy forma- 
tion model is a cosmologically consistent DM halo merger- 
tree, outlining the merging and accretion histories of DM 
halos. For the evolution of baryonic structures within these 
merging trees, we utilize the model of C09, with several 
modifications. This requires us to partition each mass ac- 
cretion history (MAH) into two phases: A 'fast accretion', 
merger dominated phase corresponding to spheroid-SMBH 
co-evolution, followed by a 'slow accretion', quiescent phase 
allowing for disk structure to form around the pre-processed 
halos (see also Zhao et al. 2003, Mo & Mao, 2004), the re- 
sulting galaxy properties at z = are thus a result of the 
sub-grid baryonic recipes, and the mass accretion histories 
of their host DM halos. Within this work we also expand 
the previous model by including several additional effects 
which are thought to be of crucial importance the low-mass 
gas-rich systems, namely, we model the effects of an ioniz- 
ing ultraviolet (UV) background, cold accretion flows and 
exploit a two-phase ISM, where star formation is governed 
by the surface density of H2 ■ 

In summary, we utilise the most up-to-date observations 
of the stellar and HI mass functions in order to derive rela- 
tionships between the host DM halo and the galaxy stellar 
and gas properties, we interpret these results using a phys- 
ically grounded SAM, and analyse the nature of the results 
by generating three model realisations, incorporating sev- 
eral contemporary recipes and highlighting how each helps 
alleviate previous model tensions. The outline of the paper 
is as follows: In §2 we present our SAM, outlining the basic 
framework from C09, and including the modifications and 
improvements, in §3 we use current observational results of 
the stellar and HI mass functions and use a numerical pro- 
cedure in order to derive relationships between host DM halo 
mass and the galaxy stellar and HI components. In §4 We 
present the results of both the observational determinations 
and the theoretical model comparisons. We conclude in §5 
by summarising the main outcomes of this work, discussing 
the implications, and the limitations. 

Throughout the paper we adopt the standard ACDM 
concordance cosmology, as constrained by WMAP 5-year 
data (Spergel et al. 2007). Specifically, we adopt a flat cos- 
mology with density parameters Qm = 0.27 and Q,a = 0.73, 
and a Hubble constant Ho = 70 km s _1 Mpc -1 . 



2 PHYSICAL MODEL 

We follow the galaxy evolution using a semi-analytical ap- 
proach, deriving analytical estimates wherever possible to 
aid simplicity. We track the evolution of dark matter, hot 
gas, cold gas and stellar mass using several parameteriza- 
tions, based on the model described in C09 but with several 
important modifications. 

Namely, we now trace the sequential infall of material 
onto both spheroid and disc structures and we include the 
effects of an ionizing UV background, adiabatic contraction 
of the DM halo and a two-phase ISM. Here we outline the 
basic framework for this model. 

We initially describe the dark matter halo evolution us- 
ing an extended Press-Schechter (EPS) algorithm based on 
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Figure 2. The derived relationships between virial mass and; 
stellar (solid/black), HI (dashed/red) and baryonic (dotted/blue) 
components. Shaded regions corresponding to the 1 — a gaussian 
errors associated with the observational determinations. 



that developed by Parkinson et al. (2008). These mergertrees 
have been tuned in order to reproduce the statistics of halo 
merger and accretion activity obtained from N-body simu- 
lations of structure formation (Springel et al. 2005). Within 
our model, we utilise these mergertrees by extracting the 
mass accretion history (MAH), i.e. the evolutionary path of 
a typical dark matter halo, obtained by tracking the most 
massive progenitor at each fragmentation event whilst mov- 
ing from z — to progressively higher redshifts (see van den 
Bosch, 2002, C09). 

Specifying the DM halo properties, we define the virial 
radius r v i r as that of a spherical volume enclosing an over- 
density in a standard way (Bryan & Norman, 1998). We 
compute the density profile for the dark matter halo using 
the fitting function of Navarro, Frenk & White (1997, here- 
after NFW): 



(1) 



In order to define the scale radius r s for our NFW profile we 
introduce the concentration parameter, which is defined to 
be c(z) = rvir/?v This quantity has been studied by several 
authors [see Bullock et al. 2001, Zhao et al. 2003 (hereafter 
Z03), Maccio et al. 2007] who found large scatter for a fixed 
halo mass, but to scale generally with the MAH of the halo. 
We adopt here the z-evolution of Z03, 



[ln(l + c) - c/(l + e) 



oc H 



z) 2a M v 



(2) 



where a is a piecewise function which can be found in Z03 
and where the normalization can be fixed using the expres- 
sion given by Maccio et al. 2007 at z = 0: 



log c = 1.071 - 0.098 



log 



12 



(3) 



From this, the scale density may be computed for a general 
profile to be p s = M vil (z) /[4nr^ /(c)], with 



/(c)=ln(l + c) 



1 + c 



(4) 



Finally, we specify the angular momentum properties 
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Figure 3. The evolution of typical model galaxies; a low-mass disk-dominated galaxy (left-hand panel), an Lq galaxy with both 
bulge and disks present (middle panel), and a high mass, spheroid dominated galaxy (right panel): DM mass (solid/black), hot 
gas (long-dash/brown), disc gas (short-dash/blue), disc stars (dotted/purple), bulge gas (triple-dot-dashed/green), bulge stars (dot- 
dashed/aquamarine), BH mass (solid/red) 



of each halo through the spin parameter, defined to be A = 
J vir i?V r 2 M^/ 2 G _1 , where i? v j r and J v i r are the total energy 
and angular momentum of the halo. Assuming that a DM 
halo acquires its angular momentum through tidal torques 
with the surrounding medium, A remains constant. It has 
been shown (Cole & Lacey, 1996) that the spin parameter 
varies little with cosmic epoch, halo mass, or environment 
and for a sample of haloes is well fitted by a log-normal 
distribution (where A = 0.04 and a — 0.5). However, for 
simplicity, here we assume that each halo has a parameter 
value of A = 0.04 which remains constant throughout the 
DM halo evolution. 

Several studies have shown the concentration evolution 
to be strongly correlated to the MAH of the halo (Z03, Li et 
al., 2006, Lu et al. 2006), finding that DM halos generally ac- 
quire their mass in two distinct phases; an initial phase char- 
acterised by rapid halo growth through major merger events, 
where the halo core structure forms, causing the gravita- 
tional potential to fluctuate rapidly, followed by a slower, 
more quiescent growth predominantly through accretion of 
material onto the outer regions of the halo. These two differ- 
ent modes are reflected in the evolution of the concentration 
parameter, which remains roughly constant during the 'fast 
accretion' phase and steadily increases during the 'slow ac- 
cretion' phase. The transition redshift zt between these two 
phases can therefore be calculated using the expression for 
the concentration parameter evolution given in Z03: 

[ln(l + Q)-c f /(l + Ct )]c- 3Q 





2a 


' M vir (z t ) 


H ° _ 




Af v ir,0 



accretion' phase, angular momentum may be readily lost by 
rapidly merging clumps of material, in an implicit merger 
scenario, resulting in the formation of a spheroidal structure, 
followed by the quiescent dissipationless infall of material in 
order to form discs (see C09 for a more detailed analysis and 
justifications of this). At difference with the prescriptions 
outlined in C09, here we expand our model to include the 
symbiotic infall of baryonic material as the DM halo evolves, 
thus 



M inf = /coiiM vir 



[ln(l + c )-co/(l + Co )] c -^ 

where H{z) is the Hubble radius, '0' denotes quantities eval- 
uated at the present cosmic time, and c t = 4 and a = 0.48 
are scaled to match N-body simulations as in Z03. 

By associating the two phases of DM evolution to two 
growth mechanisms for the baryonic sector, the fast and 
slow phases give rise to the formation of bulges and discs 
respectively (Mo & Mao, 2004, C09). Since within the 'fast 



,(5) 



(6) 



We include the effects of an ionizing radiation background 
taking the prescriptions outlined in (Gnedin et al., 2004, 
Somerville et al., 2008), which is able to partially reduce the 
baryonic content in low-mass systems, thus 

rn 



/coll(-Mvi r 



(l + O^Af/^/Afvir) 3 ' 



(7) 



where Mj(£) is the filtering mass at a given redshift, com- 
puted using the equations in (Kravtsov, Gnedin & Kyplin, 
2004, Appendix B). A second improvement over our previ- 
ous model (C09) is to include the effects of cold accretion 
flows, shown to be the predominant mechanism leading to 
the formation of low-mass systems. Below a critical mass 

M c = max[l, lO 1 ' 3 ^^] , (8) 

where M s = 2 x 1O 12 M0 and z c = 3.2, we assume that all 
gas accreted onto DM halos is not shock heated to the virial 
temperature of the DM halo, but streams in on a dynamical 
time (see Dekel et al., 2008, 2009, Cattaneo et al, 2006). We 
note that below the shock heating mass scale it has been 
shown that rapid cooling does not allow for the formation 
of a stable virial shock (Keres et al., 2005) resulting in gas 
flowing unperturbed into the central regions of the DM halo. 
Thus, in halos below this mass the collapse happens on the 
dynamical timescale of the system (t co il = *dyn), whereas in 
halos above this mass t co l\ = max[£dyn, icooi], where the cool- 
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Figure 4. Stellar mass function output by each of our three model 
realisations compared with the Bell et al. 2003a estimate. We find 
that the most sophisticated treatment of the baryonic physics re- 
sults in agreement to the observations in all but the lowest mass 
halos (M a < ICPMq), where we find that we overproduce the 
stellar matter by approximately a factor two. Note that, due to 
large uncertanties in observational mass- to-light ratios, we asso- 
ciate 0.3dex errors in the stellar mass determinations, and 0.4dex 
in the value of <p(M B tars)i shown by the shaded area. Model error- 
bars represent poissonian uncertanties due to the finite sample 
size of synthetic galaxies. 



ing timescale t coo \ is computed in a standard way, assuming 
material is shock heated to the virial temperature. The ef- 
fects of this cold accretion is to moderately enhance star 
formation at high redshifts relative to the scenario where all 
material is shock heated, in this way, we model the infall 
and cooling dynamics using the most current recipes (see 
Somerville et al, 2008). 

In order to calculate the cooling time of the hot gas 
phase, we assume an isothermal gas in hydrostatic equilib- 
rium within the NFW profile, such that 



Phot( 



Po exp 



-27 



fill 



ln(l + r/r 3 ) 
r/r s 



with 



(3 = 



9>-Kp,m v Gp s r a 
27k B T viT 



(9) 



(10) 



where T v i r is the virial temperature, p is the mean molec- 
ular mass and po is calculated by normalizing to the total 
hot gas mass at any given redshift. Also, we assume that the 
cooling function A(T, Z) is given by the tabulated function 
of Sutherland & Dopita, 1993, and assume that the infalling 
baryonic material is unprocessed and therefore has primor- 
dial metallicity, Z = 1O~ 3 Z . 



the overall dynamical time of the forming halo. Therefore, 
within the 'fast accretion' phase we neglect the effects of 
angular momentum of cooling proto-galactic gas, which will 
result in the formation of a spheroidal gaseous system at a 
rate 



Afcoii(z) =4tt 



I "vir( z ) 



r- 2 phot(?-, z) 
tcou(r, z) 



dr , 



(11) 



where t co \\(r, z) is determined by the cold-accretion recipe 
of Eqs. 8, 9 and 10. We denote this cold gas spheroidal com- 
ponent ('gaseous bulge'), which act as a reservoir for star 
formation, by Mb tgas (t), and stress that besides dissipative 
collapse, it can also grow through merger events and disc in- 
stabilities. Also, our model includes other two components: a 
spheroidal stellar component ('stellar bulge') Mb tStar (t) and 
a low angular momentum cold gas reservoir M re3 (t), which 
acts as a source of material eligible to accrete onto a central 
black hole. 

For simplicitjQ, we assume that the reservoir can be 
described by an exponential disc surface density 

-r/r res (z) 



(12) 



with the scale radius r ICS being proportional to the influence 
radius of the SMBH (r rcs = aGM SMBH /l/ 2 ir , with a « 100). 
Also, since we have no a priori information about the geo- 
metric distribution of baryonic matter within the bulge sys- 
tem, and since the dynamics and thus evolution of disc struc- 
ture is correlated to the mass and geometry of the bulge 
structure, we assume that the bulge stellar and gaseous 
masses settle into a Hernquist density profile 



Pb ( 



MZ 



''!> 



27r r(r + rb) 3 



* = stars, gas , 



(13) 



where the scale radius of this profile is related to the half 
light radius by r& = 1.8152J2 c ff - Using the fitting of Shen 
et al. (2003), we take the parametrization as a function of 
bulge mass to be 



log(-Refl) = 



-5.54 + 0.56 log (M 6 ) 
-1.21 +0.141og(M 6 ) 



log(Mi,) > 10.3] 
log(Mfc) < 10.3] 



The star formation rate (SFR) per annulus in the 
gaseous bulge may be computed as 



dr 



(r, t) = Anr 



2 pb,gas(r) 

tgas(r) 



(14) 



where t gas (r) is the dynamical time for the gas in the bulge. 
Therefore in order to compute the total SFR we must inte- 
grate this expression over all radii. 

Energetic feedback due to supernova events may trans- 
fer significant energy into the cold ISM, causing it to be 
re-heated and ejected from the system. Therefore, by con- 
sidering energy balance in the ISM, supernovae feedback is 
able to remove gas from the bulge at a rate: 



2.1 Dissipative gas collapse 

During the high redshift domain, it has been shown that 
the majority of the angular momentum of the collapsing 
proto-galactic gas dissipates as it collapses and condenses to 
the centres of DM haloes (see Navarro & Steinmetz, 2000). 
Also, within C09 we show that during the early collapse 
phase, DM halo subunits merge on timescales shorter than 



2 We need to make an assumption about the reservoir geometry 
because that is needed to calculate the velocity of the composite 
system V c [needed e.g. in Eqs. 27 and 32, the adiabatic halo con- 
traction factor r [Eqs. 28 and 29] and the gravitational potential 
tf> of the composite system appearing in Eqs. 15 and 34. However, 
the geometry of the reservoir is not expected to have a major 
impact on our results, given its small size relative to the other 
components. 



6 M. Cook et al. 



Table 1. Values of the free parameters of our model. We note that the value of k aC c, despite being a free parameter, has a minor impact 
on our results 



Description 


Symbol 


Fiducial value 


Reference in the text 


Impact on this work 


SN feedback efficiency (bulge) 


<=SN,b 


0.5 


Eq. 15 


Strong 


SN feedback efficiency (disc) 


e SN,d 


0.7 


Eq. 35 


Strong 


Reservoir growth rate 


A 

sires 


W- 3 M e yr- 1 


Eq. 18 


Strong 


QSO feedback efficiency 


fh 


10" 4 


Eq. 22, 23 


Strong 


Viscous accretion rate 


kacc 


10~ 2 


Eq. 20 


Weak 


Note. - A Romano IMF (/>(m*) is 


adopted: <j. 


>(m*) oc m^" 1 ' 25 


for m* Mq and (^(m*) 


oc m~ 0A for m* ^ Mq. 



M b % s (t) = 



esN,bEsNT]SNdijjb (r, t) / dr 
4>(r,t) 



dr , 



(15) 



where t/sn is the number of Type II supernovae expected 
per solar mass of stars formeqj, Esn is the kinetic energy 
released per supernova event, and e$N,b is the efficiency of 
supernovae energy transfer used to remove the cold gas. Fi- 
nally, <f>(r, t) is the gravitational potential of the compos- 
ite system (bulge, reservoir, disc and DM). Therefore, using 
this prescription, we see that supernovae feedback is partic- 
ularly efficient in smaller halos with shallower gravitational 
potential wells but relatively inefficient in larger halos. Us- 
ing scaling relations, Granato et al. (2000) were able to show 
that stars form faster in larger systems, thus exhibiting the 
observed 'anti-hierarchical' behavior of spheroid galaxies. 

A growing body of evidence is now showing that the 
evolution of both the SFR within spheroids and the fueling 
of SMBH's are proportional to one another (Haiman, Ciotti 
& Ostriker, 2004). A proposed mechanism to account for 
this phenomenon has been discussed (Kawakatu, Umemura 
& Mori, 2003): Radiation drag due to stellar radiation may 
result in the loss of angular momentum at a rate well ap- 
proximated within a clumpy ISM by: 



dln(J) 
dt 



--'sph 



c 2 M b ,, 



"(I 



(16) 



where L sp h is the total stellar luminosity and r is the ef- 
fective optical depth of the spheroid (r = fNi nt where f is 
the average optical depth and Ni„t is the average number 
of clouds intersecting a light path). Upon loss of angular 
momentum this gas may flow into the nuclear region, gener- 
ating the reservoir of low-J material which fuels BH growth 
at the rate (Granato et al., 2004) 



1.2 x lQ- 3 4j b (t)(l - e' T )M Q yr" 1 



(17) 



We note that, within this work we assume that r is constant 
for simplicity, which allows one to rewrite Eq. 17 as 



Mr 



A res ip b (t) . 



(18) 



where A res is a free parameter. The cold gas stored in the 
reservoir is then expected to accrete onto the central SMBH 
with an accretion rate 



M bh = ram\M vi3a ,M e , 



(19) 



In this formula, the viscous accretion rate is given by 
(Granato et al., 2004) 

3 A Romano IMF <p(m t ) is adopted: 0(m*) oc m+ 1,25 for m* 
Mq and c/>(m*) oc m^ 11,4 for m* ^ Mq. This gives t)sn = 5 X 



10" 3 M, 



where k acc £= 
simply M edd 



ct 3 / M res \ 
G \ M bh ) 



(20) 



10 , whilst the Eddington accretion rate is 



L e dd/V C i wit 



10 



0.15 and 



L Edd w 1.26 x 10 



46 Mbh (i) 
10 s Mq 



erg s 



(21) 



QSO activity affects the interstellar medium of the 
host galaxy and also the surrounding intergalactic medium 
through both radiative heating and the kinetic energy input 
through gas outflows. Assuming that a fraction fh (which 
we treat as a free parameter) of the SMBH luminosity Lh is 
transferred into the cold and hot gas phases, it is possible to 
compute the amount of cold and hot gas which is removed 
from the hot gas and gasoeus bulge phases as in Granato et 
al., 2004: 



MQSO 

b,gas 



Jh- 



Mi 



b,gas 



M 



QSO 



3 (J 2 M hot + M b , gas 
2 Lh Mhot 



(22) 
(23) 



3 a 2 Mhot + M b , gas 

This material is assumed to be ejected 



where a = 0.65U v i r 
from the system. 

For the chemical evolution of the cold bulge gas, we use 
the simple instant-recycling approximation (IRA), whereby 
a fraction of mass is instantly returned into the cold gas 
phase in the form of processed materia^]. In particular, this 
implies that the effective SFR which enters the evolution 
equations for the gas and star bulge masses is given by 



Mi FR (t) 



(1-R) 



(r,t)dr , 



(24) 







with 7? = 0.25. Also, we assume that Mhydrogen = 
0.71M co i d , where the factor takes into account the contri- 
bution of Helium and other heavier elements. 



2.2 Dissipationless gas collapse 

By assuming that material may collapse to form a cool 
gaseous disc structure during the slow accretion phase, we 



4 This value of r] corresponds to rapidly spinning SMBH with 
spin parameter a Kb 0.9 (Bardeen, 1970). 

5 We adopt a Romano et al., 2005 IMF, which has the standard 
Salpeter slope 1.25 in the high mass tail, and flattens to a slope 
0.4 below 1 Mq. As shown in Romano et al. (2005), this performs 
better than the Salpeter one in reproducing the detailed chemical 
properties of elliptical galaxies. 
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Figure 5. Baryonic mass function output by our three model 
realisations as compared to the Bell et al. 2003b estimates. As 
is shown, the 'standard' approach over predicts significantly the 
amount of baryonic material at all masses, whereas adding the 
effects of a UV background and then a two-phase ISM results 
in progressive improvements. The most sophisticated realisation 
results in a good fit to the observational mass function across all 
observable masses. Note that the shaded area represents the 1 — a 
errors, discussed in fig. 4. Model error-bars represent poissonian 
uncertanties due to the finite sample size of synthetic galaxies. 



may add material to the disc structure at a rate which is 
given, as in the spheroidal case, by 



Mcoii(z) =4tt 



r(2> r 2 p hot (r, z) 



tcoi\(r, z) 



dr , 



(25) 



where t co u(r,z) is determined again by the cold-accretion 
determinations (see Eqs. 8, 9 and 10). If we assume a dissi- 
pationless collapse of material upon cooling within the dark 
matter halo, we may relate the dark halo virial radius and 
spin parameter to the forming disc scale radius. In particu- 
lar, if we assume an exponential disc surface density profile 
for the stellar and gaseous components, 



2%{r,z) = Y,* (z)e 



-r/ r d(z) 



* = stars, gas , 



(26) 



the disc scale radius r d {z) evolves according to the scaling 
r d (z) = (2T:y 1 ' 2 (j d /m d )Xr vir (z)f(cy 1 ' 2 f r (X,c,m d j d ). The 
function /(A, c,m d j d ) may be exactly determined through 
(Mo, Mao & White, 1998) 



/(A,c, m d ,j d ) 



i 2 V c (r d u) 
u tt"t fdu 



V c (r vil 



(27) 



where V c (r) is the velocity profile of the composite system 
(bulge, reservoir, disc and DM) and where m d and j d are 
the ratios between the total mass and angular momentum 
of the disc component and the DM halo mass. More specif- 
ically, we take m d = (M,f ars + Mf as )/M vir , and we assume 
j d = m d (Mo, Mao & White, 1998). In order to account for 
adiabatic halo response, we take the standard prescription 
of Blumenthal (1986). In particular, denoting by Mx(r) the 
mass of the a given component 'A" enclosed by a radius r, 
from the angular momentum conservation one obtains 



Mi(n)ri = Mf(r f )r f 



(28) 



where r, and rj are respectively the initial and final radius of 
the shell under consideration, the initial mass distribution 
Mi(ri) is simply given by the NFW density profile, while 
Mf(rf) is the final mass distribution. Also, mass conserva- 
tion easily gives 

M f (r f ) = M d (r s ) + M b (r f ) + M DM (r f ) + M les (r f ) = 
M d (r f ) + M b (r f ) + M las (r f ) + (1 - f gal )Mi(n) , (29) 

where f ga i = M ga i/M viT (with M gat = M d + M b + M TCE ). By 
assuming spherical collapse without shell crossing, one can 
adopt the ansatz ry = Tr^, with T = const (Blumenthal, 
1986), and Eqs. 28 and 29 may be solved numerically for 
the contraction factor T. 

When the surface density of the gaseous disc increases, 
the cold gas becomes available to form stars. However, at 
present, star formation is poorly understood from both a 
microscopic, and large-scales. Therefore we parameterise the 
star formation using an empirical Schmidt law (Kennicutt, 
1998) whereby the star formation rate is related to the sur- 
face density of cold disc gas: 



E sfr (r, z) = e„f 



vigas 



(r, z) 



M Q pc- 



M Q kpc 2 y r 1 



(30) 



where e B f = 2.5 x 10~ 4 controls the star formation efficiency 
and n = 1.4 is fixed to match the properties of spiral galax- 
ies. This law is empirically proven to hold over many orders 
of magnitude of gas surface density, but was shown to break 
down at large radii (Kennicutt, 1998) where the surface den- 
sity drops below a critical value, roughly corresponding to 
the Toomre (1964) stability criterion. We take this critical 
surface density to be 



E c (r) = 



3.36GQ' 



(31) 



Toomre, 1964. Where a g = 6km s~ is the velocity disper- 
sion of the gas, Q = 1.5 is a dimensionless constant and n(r) 
is the epicyclic frequency, given by 



dV c (r) 



V c (r) dr 



1/2 



(32) 



Therefore, the conversion rate of gas mass to stellar mass is 
computed as 



i/j d {z) = 2-KEst / r£2, gas (r, z)dr, 



(33) 



where r c may be calculated by solving E| as (r c ,z) = S c (r c ) 
for r c . 

In order to account for feedback due to supernovae 
events, we may compute the amount of cold gas which is 
ejected from the system at each disc radius. In order to 
remove this cold gas from the disc, the supernovae feed- 
back must be sufficient to unbind it, therefore we compare 
the amount of energy released through supernovae events at 
each disc radius with the binding energy at the same radius: 



Ssn(t, z) 



£SN,d-EsN?;sivS s f r (r, z) 
4>(r,z) 



(34) 



where </>(r, z) again, is the binding energy of the composite 
system (bulge, disc, reservoir and DM). The total amount 
of cold gas ejected from the system is then given by 
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Figure 6. The HI mass function produced by our model as com- 
pared to the Zwaan et al. 2005 HIPASS galaxy sample. Using the 
combined effects of a UV background photo-ionizing radiation 
and a two-phase ISM we suppress the formation of low mass sys- 
tems. We find that we reproduce both the normalisation and high 
mass cutoff accurately, but still slightly overproduce the number 
systems with Afjjl < 10 9 Mq, a minor effect discussed in §4 & 5. 
Observational error bars are associated to 1 — a uncertanties in de- 
termination of HI masses. Model error-bars represent poissonian 
uncertanties due to the finite sample size of synthetic galaxies. 



Mg N (z) = 2ir 



rEsN(r, z)dr 



(35) 



For the chemical evolution of the cold disc gas, we use 
again the IRA approximation, and the 'effective' SFR which 
enter the evolution equations for the gas and star disc masses 
is given by 

d 



SFR 



(z) = (1 - R)Mz) ■ 



(36) 



with R — 0.25. Also, we assume again that Mhydrogen = 
0.7lM co i d . 

Finally, it is known that when discs become self- 
gravitating they are likely to develop bar instabilities, get 
disrupted and transfer material to the spheroidal compo- 
nent (Christodoulou, Shlosman & Tohline, 1995). We there- 
fore assume that a stellar or gaseous disk is stable if 



stars , gas , 



(37) 



K(2.2r d ) 
(GM* isc /r d )V2 

where atgg* = 1.1 and agf t = 0.9 [see (Mo, Mao & White, 
1998) and references therein]. If we find that discs become 
unstable, we assume they get disrupted in a dynamical time 
and transfer their material (either stars or gas) to the bulge 
components. 



2.3 Improved star formation law 

Assuming that star formation may only take place inside 
dense molecular clouds several authors have shown that the 
Schmidt-Kennicutt star formation law (Eq. 30) may be re- 
produced within large mass systems by assuming that the 
star formation rate is proportional to the molecular cloud 
mass (Blitz & Rosolowski, 2006, Dutton & van den Bosch, 
2009), thus: 



Where e sf = 13Gyr _1 and £ mo i,HCN = /moi-RHCN is the 
molecular mass surface density as traced by HCN (see Gao 
& Solomon, 2004, Wu et al. 2005). Thus, calculating the 
ratio of molecular gas to atomic gas allows us to compute 
the star formation rate at all scales. The fraction of gas in 
discs which is molecular has been extensively analysed, and 
shown to be closely related to the mid-plane pressure within 
discs (Blitz & Rosolowsky, 2006), given by: 



-GEg (Eg + (<7g/<7 3 )E s 



(39) 



Where, following the detailed prescriptions of Dutton & 
van den Bosch, 2009, assuming a constant cr g /a s — 0.1. Re- 
lating the mid-plane pressure to the formation of molecular 
clouds yields: 



Hn 



4.3 x 10 4 



(40) 



Thus, the molecular fraction is given by / mo i = 
-Rmoi/(-Rmoi + 1), in order to relate this to the HCN frac- 
tion, we must further compute the fitting relation of Blitz 
& Rosolowski, 2006: 



CN 



0.1*(l + £ /(2OOM pc 



- 2)) 0.4 



(41) 



Therefore, we find that, in the high mass (and thus den- 
sity) galaxies, where the molecular fraction is typically ~ 1, 
we recover the standard Schmidt-Kennicutt star formation 
power law, with exponent 1.4, whereas in the low-density 
galaxies we asymptote towards an exponent of 2.84, sup- 
pressing star formation in these systems, in accordance with 
observations (see Dutton & van den Bosch, 2009 and Blitz & 
Rosolowski, 2006 for a detailed description) . Thus, through- 
out the evolution of each galaxy, we partition the ISM into 
HI and H2 components using the above relations, and com- 
pute the star formation law (and therefore supernovae feed- 
back) using equation 42 in order to self-consistently model 
each galaxy under this improved star formation law. 



3 THE HI MASS RELATIONSHIPS 

In order to investigate the relationship between the stellar 
and the gas component in late-type galaxies, we follow the 
procedure already exploited by Shankar et al. (2006), and 
discussed in detail in Evoli et al., 2009. We defer the reader 
to these papers for a more detailed discussion and highlight 
the main results here. If two galaxy properties, q and p, obey 
a one-to-one relationship, we can write 



dp 

— dq = tp(q)dq 



(42) 



where ^(q) is the number density of galaxies with mea- 
sured property between q and q + dq and 4>(p) is the corre- 
sponding number density for the variable p. The solution is 
based on a numerical scheme that imposes that the number 
of galaxies with q above a certain value q must be equal to 
the number of galaxies with p above p, i.e., 



(43) 



4>(p)dp= I tp(q)dq 
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Figure 7. The H2 mass function produced by our model as com- 
pared to the estimations by Keres, 2003. We find that we re- 
produce both the normalisation, low mass, and high mass cutoff 
within the observational range. Observational error bars are as- 
sociated to 1 — a uncertanties in determination of HI masses. 
Model error-bars represent poissonian uncertanties due to the fi- 
nite sample size of synthetic galaxies. 



In the following the HI mass function (HIMF) is given 
by; p = Mm and (j>{p) = HIMF(Mhi), while the variable q is 
the stellar mass (M*) and ij)(q) the corresponding Galactic 
Stellar mass function (GSMF(M*)). 

Using the HI Parkes All-Sky Survey HIPASS (Meyer et 
al. 2004) it has become possible to map the distribution of 
HI in the nearby Universe. The HIMF has been fitted with 
a Schecter function with a power index of a = — 1.37T0.03, 
a characteristic mass oflog(MH J /M s ) = 9.80±0.03/if 5 2 and 
a normalization of 9* = (6.0±0.8) x W~ s hj 5 Mpc" 3 dex" 1 . 
Nakamura et al. (2003) estimated the LF in the r* band 
for early- and late- type galaxies separately. Shankar et al. 
(2006) used the results of Zucca et al. (1997) and Love- 
day (1998) to extend these results to lower luminosities and 
giving a good fit for the late-type galaxy LF in the range 
10 7 L Q < L r < 3 • 1O U L . The GSMF, holding over the 
mass range 10 s < M* < 10 12 is compared with the one of 
Bell (2003a), and is in good agreement within the uncertain- 
ties due to the MLR, estimated to be around 30%. 

The correlation of galaxy properties with the halo mass 
(Mh) is extremely relevant in the framework of galaxy for- 
mation theories. To constrain such a relation we used the 
association between the stellar mass and the host halo mass 
derived in Shankar et al., 2006 (Eqn.12 & Fig.l therein), 
which has been obtained with the method already mentioned 
but using either for the stellar mass function and the galactic 
halo mass function a fitting of the observations/simulation 
for all the galaxies. Finding relations between each of the 
mas components within galaxies is of great importance, be- 
cause it allows for constraints on the direct outputs from 
models, not requiring significant post processing, and pro- 
viding useful links between the dark and luminous compo- 
nents of galaxies. We find that the results of the numerical 
outputs may be well approximated by the following analytic 
fitting functions: 



Mhi 



(M vir /5.7 x 1O 1O M ) 5 



(44) 



(45) 



6.07 x 1O 8 M 1 + (M vir /5.7 x 1O 1O M ) 4 - 76 

Fig. 2 shows the ratio of the baryonic mass components 
(stellar, gas and total) to the initial baryon mass associ- 
ated with each DM halo obtained using the above-mentioned 
methods (and explained in Evoli et al., 2009), illustrating the 
inefficiency of galaxies, especially those of low halo mass, in 
retaining baryons. These relationships can be used to tightly 
constrain theoretical galaxy formation models in order to 
interpret the physical processes relevant to shape these re- 
lations. 



4 RESULTS 

Throughout this section we present the results of model 
realisations, each comprising of a sample of w 10 3 galax- 
ies in logarithmic virial mass increments in the mass range 
9.5 < log(M v (z = O)/A'f ) < 13.5 in order to encompass the 
observational range of galaxy-hosting DM halos. In order 
to make statistical predictions with each generated galaxy 
sample at z = we weight each DM halo with the GHMF, 
which mirrors the Sheth & Tormen (2002) mass function 
within most galaxy sized DM haloes, but is derived in order 
to account for the increasing probability of multiple galaxy 
occupation in the highest mass haloes (see Shankar et al. 
2006). Throughout this section we also plot the results for 
'early' and 'late' type galaxies, by parting the populations 
into [M bulge /M tota {\ > 0.5 and [M bulg( ,/ M tota {\ < 0.5 respec- 
tively. 



4.1 Component evolution 

In order to illustrate the general behavior of our galaxy for- 
mation framework, Fig. 3 highlights the evolution of each 
mass component from 2 = 8 up to the present epoch. In 
the left-hand, center and right-hand panels we show a low- 
mass, intermediate mass and high-mass galaxy respectively, 
in order to highlight typical model outputs for each system. 
As can be seen, the evolution of each galaxy differs signifi- 
cantly due both to scatter at each mass through monte carlo 
selected mass accretion histories, and to the relative differ- 
ences in efficiencies of the competing processes of infall, star 
formation and feedback on different mass-scales. 

We find, in broad agreement with observations, that 
a typical low-mass galaxy (with M v (z = 0) « 10 10 M ) 
supports the growth of a disc structure from high redshifts 
showing extended star formation up to z = resulting in a 
gas-rich disc dominated galaxy with a negligible bulge com- 
ponent. In intermediate mass, L galaxies, at high-redshifts 
gaseous collapse onto a spheroid system results in the co- 
evolution of SMBH's and the spheroidal component result- 
ing in a gas poor, 'red and dead' spheroidal stellar compo- 
nent which acts as the bulge within the resultant formed 
galaxy as the disc component grows steadily from z ^ 1.5. 
The high mass galaxy shows rapid early growth of the SMBH 
and spheroid component, followed by an epoch of dormancy 
in the star formation, no significant disc structure may form 
in these systems and ongoing star formation is prevented due 
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Figure 8. Deriving the relationships in §3 we show the baryon 
budget in the form of stellar material in DM halos of different 
mass. Only by exploiting the effects of a UV background and a 
two-phase ISM can we reproduce the global form of this rela- 
tionship, finding good agreements to observations at all masses, 
with a slight over-prediction only in the lowest mass halos. As 
can be seen, 'standard' approaches significantly overproduce stel- 
lar material on all mass ranges, but the discrepancy is clearest 
in the low mass systems. The shaded region represents the f — a 
observational uncertanty. Model error-bars represent poissonian 
uncertanties due to the finite sample size of synthetic galaxies. 



to the presence of a large SMBH which acts to effectively 
expel any residual gas which may infall later, only allow- 
ing for extremely mild star formation following the main 
growth phase. The resulting galaxy spheroidal component 
thus comprises an old stellar population, with little gas and 
negligible star formation, and the disc component is over 
an order of magnitude smaller than the bulge, and gas rich, 
showing little star formation. 

These results are consistent with the evolutionary histo- 
ries tuned to match the chemical properties of local galax- 
ies of different morphological type and different host halo 
masses (Calura, Pipino & Matteucci, 2008 & references 
therein) which have also been shown to be consistent with 
photometry within the local Universe (Schurer et al. 2009) 
and are generally in agreement with observations of statis- 
tical samples of galaxies (Driver, et al. 2006), however, as 
can also be seen, due to the stochastic nature of the model, 
fluctuations in galaxy properties are expected, and thus we 
must revert to statistically representative samples of galaxy 
populations in order to make more robust comparisons. 

Within these systems in-falling baryonic material is 
shock heated and therefore may form a static atmosphere 
of hot-gas. Feedback from a growing SMBH may thus halt 
the cooling of this gas, quenching star formation and sweep- 
ing out the ISM (See Granato et al., 2004). Within these 
systems it is thought that baryonic material is shock heated 
and therefore cooling and feedback from a growing SMBH 
may quench star formation by sweeping out the ISM (see 
Granato et al. 2004) 



4.2 Stellar and Baryon mass functions 

The local stellar and baryonic mass functions provide a pow- 
erful constraint on theoretical models of galaxy formation: 
Encompassing much of the relevant physical processes which 
determine the assembly of baryons within DM halos. 

In order for models to reproduce observational results, 
it has become clear that physical processes of gas accretion, 
supernovae feedback and photo-ionization are most effec- 
tive in the lowest mass DM halos where the shallow poten- 
tial wells are inefficient at trapping and holding baryonic 
material. Thus these processes drive the evolution of the 
faint-end slope of the mass functions (Benson et al., 2002), 
whereas the brightest galaxies (above L* ) are embedded 
within large DM halos which effectively trap baryonic ma- 
terial. Within these systems in-falling baryonic material is 
shock heated and therefore may form a static atmosphere of 
hot-gas. Feedback from a growing SMBH may thus halt the 
cooling of this gas, quenching star formation and sweeping 
out the ISM (See Granato et al., 2004). Coupled with the 
increasing subhalo contribution (through the cluster mass 
function) whereby in DM halos with M v (z = 0) > 10 13 the 
probability of a single galaxy occupation is low, gives rise to 
the relatively sharp cutoff of stellar mass within the largest 
DM halos (see Shankar et al. 2006 & Somerville et al. 2008, 
& contained references). 

Rather than computing the spectral energy distribu- 
tion (SED) assuming any mass-to-light ratio for each model 
galaxy (relying upon further model assumptions such as 
the dust-to-gas ratio, molecular cloud structure and opti- 
cal depth etc) by comparing the stellar and total baryon 
budgets within each DM halo, we provide the most direct 
analysis of model outputs. However, we also note that un- 
certainties in the observational conversion of luminosity to 
stellar (and total baryonic) mass are systematically related 
to the spectral energy SED fitting methods used in order to 
extract physical quantities from the multi-wavelength obser- 
vations, and on the quality of the observations themselves, 
this is known to have large uncertanties and we hope in a 
subsequent work, to utilize synthetic spectra using the de- 
tailed star formation histories, galaxy geometries, gas and 
dust content, in order to self-consistently model the multi- 
wavelength SED and make comparisons with the luminous 
properties of galaxies. In Fig. 4 we show the Schecter func- 
tion fit to the Bell et al. (2003a) estimate for the local stellar 
mass function, in order to generate this mass function the 
authors utilize a large sample of galaxies from the Two Mi- 
cron All Sky Survey (2MASS) and the Sloan Digital Sky Sur- 
vey (SDSS) converting galaxy luminosity into stellar mass 
using simple models to convert the optical and near infrared 
observations into stellar masses. 

As we have shown in Fig. 4, using three levels of so- 
phistication we are able to highlight the differences be- 
tween; a 'standard' model, whereby we ignore the effects of a 
photo-ionizing background and employ a standard Schmidt- 
Kennicutt (1998) star formation law, a 'UV model whereby 
we include the suppression of mass flowing into galaxies 
through photo-ionization (see Eqn. 6 & 7), and a 'sophisti- 
cated' model, which combines the effects of a photo-ionizing 
background with a two-phase ISM and modified star forma- 
tion law. We will refer to these model names throughout the 
next sections. 
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More specifically, using a 'standard' model we signifi- 
cantly over-produce the number density of low mass galax- 
ies. This occurs even when we increase supernovae feedback 
efficiencies to extremely high levels (see Mo et al. 2005 for a 
discussion), further pushing the parameter to higher values 
than e sn — 0.7 would result in largely un-physical models, 
and it is clear that reproducing observations using this basic 
framework does not occur naturalhjf]. Secondly, increasing 
the supernovae feedback efficiencies whilst lowering the star 
formation efficiencies simply has the effect of reducing the 
stellar mass but not the gaseous mass, resulting in gas frac- 
tions which increase with increasing supernovae efficiency 
and thus are in disagreement with observations (see Mo et 
al., 2005 for a detailed discussion). 

By suppressing the initial in-fall of material due to an 
ionizing UV background radiation, we are able to improve 
agreements with the mass function. However, we achieve the 
best reproduction of the stellar mass function through the 
additional reduction in star formation efficiency when em- 
ploying the 'sophisticated' star formation law, which is de- 
temined by the amount of H2 gas present in galaxy discs. 

Alone, the stellar mass function provides a crucial ob- 
servation for any physical galaxy formation model to re- 
produce, but by simultaneously comparing both the stel- 
lar and gaseous properties of galaxies, we are able to break 
degeneracies between gas infall and cooling, star formation 
and feedback processes, resulting in a significantly improved 
constraint on the theoretical framework. This has led pre- 
vious authors to claim that the standard model, whereby 
low mass systems support efficient gas cooling leading to 
large gaseous rotationally supported discs, which then un- 
dergo mild star formation, significantly over-predicts galaxy 
HI masses when compared to observations (see Mo et al. 
2005) causing serious tension between theory and observa- 
tion. 

In Fig. 5 we compare model results with the Bell et al. 
(2003b) baryon mass function. We find that suppressing the 
infall of baryonic gas due to an ionizing UV background 
significantly improves agreement between observation and 
model output, since material is prevented from in-falling 
into the halos initially and thus requiring less feedback in 
order to gain agreement of the low stellar and gas fractions 
in these halos. A final, and further improvement between 
model and observation is achieved by partitioning the ISM 
into neutral and molecular gas, which has the effect of reduc- 
ing the SFR efficiency preferentially in the lowest mass DM 
halos. We find that the best agreement between model and 
observation naturally results from using the most sophisti- 
cated treatment of the ISM physics, and the initial infall 
of gas. We do however, still find discrepancies in the lowest 
mass halos (M st ar(z = 0) < 10 9 ), over-predicting the baryon 
and stellar content within these DM halos, and slightly over- 
producing the number density of the most massive galaxies. 
This issue is discussed in §5. 

6 We note also that we use a universal stellar initial mass func- 
tion (IMF) which has been constrained in order to reproduce 
the chemical properties of galaxies accurately, increased super- 
novae rates are achievable using a more top-heavy IMF, however 
this has a strong chance to offset other properties of the formed 
galaxies with respect to observations (see Romano et al. 2005 for 
a discussion.) 
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Figure 9. The fraction of baryonic matter in the form of HI 
relative to the cosmological fraction showing that using the 'so- 
phisticated' model we are able to successfully reproduce the de- 
rived relationships in §3 accurately down to relatively low halo 
masses, but over-predict the HI mass by a factor of two within 
the lowest mass halos, we discuss the implications of this in §5. 
The shaded region represents the 1 — a observational uncertanty. 
Model error-bars represent poissonian uncertanties due to the fi- 
nite sample size of synthetic galaxies. 

However, the observations (with their large uncertain- 
ties), remain relatively insensitive to the detailed properties 
of the galaxies at different masses, a strong motivation for 
the method adopted in §3, whereby linear scales amplify the 
discrepancies between models. 

4.3 HI and H2 mass functions 

Due to the implemented star formation law within the 'so- 
phisticated' model, whereby we partition the total cold gas 
mass into atomic and molecular HydrogerQ, where f mo i — 
Rmoi/(R mol +i) and R mo i is given by Eq.40 at each time- 
increment in order to self-consistently generate a two phase 
ISM, we are able to make simultaneous predictions for both 
the HI and H2 mass functions. In Fig.6 we compare our 
model results to the observational estimates of Zwaan et al. 
2005. The sample, comprising of 4315 extragalactic emis- 
sion line estimated HI masses from the HIPASS catalogue, 
is complete down to Mm ~ 10 8 Mq, which allows for an un- 
precedented determination of the low mass slope. We find 
that above Mhi ~ 10 9 Mq we accurately reproduce the form 
and normalisation of the mass function, but we do signifi- 
cantly over-produce the HI masses at the low-mass extreme. 

Secondly, in Fig. 7 we compare the H2 mass function as 
derived from the CO-mass function determination by Keres 
et al. 2003 based on the FCRAO extragalactic CO survey of 
200 galaxies (Young et al. 1995). The conversion from CO to 
H2 is achieved through a '^-factor' which proves to be a del- 
icate and difficult task to determine, we utilize the H2 mass 
function derived by Obreschkow & Rawlins, 2009 as this 
appears to be the most robust determination, and also plot 

7 Taking Mh 2 = fmolM hydrogen and M HI = (1 - 

fmol)M hydrog 

en 



12 M. Cook et al. 



the H2 mass function derived within the originalKeres et al., 
2005 work (using a constant x-factor). We refer the reader 
to their work for further details. We find that within obser- 
vational errors, we reproduce this function over the entire 
observational range, however, we note that the large uncer- 
tanties within the determination of the precise value of the 
X-factor means that we are relatively loosely constrained, 
and we view this result as a general prediction of our model, 
rather than a constrained observational match. 

Physically interpreting these results, we may conclude 
that utilising a two-phase ISM whereby star formation scales 
only indirectly with the total gas mass, but directly with the 
molecular H2 mass, we find that we may accurately repro- 
duce both the HI and H2 mass functions of galaxies apart 
from an over-prediction of neutral HI within the lowest mass 
halos. This discrepancy also appears within the stellar mass 
function (Fig.4) and the baryonic mass function (Fig. 5), al- 
though it is greatly improved using the 'sophisticated' ISM 
treatment. 

4.4 Stellar-to-Dark matter properties 

Exploiting the methods outlined in §3 we are able to make 
observational mappings between DM halo masses and their 
contained stellar mass. In Fig. 8 we show the relative effi- 
ciency of conversion of baryonic material into stellar mass 
within DM halos relative to the cosmological baryon to 
dark matter ratio ft = — 0.16. These relation- 

ships provide direct observational constraints onto galaxy 
formation models and allow for a direct analysis of the rel- 
evant processes occurring over a range of mass scales (see 
Shankar et al. 2006). As can be seen, the lowest mass systems 
(log(M v /Mo) < 11) are strongly dark matter dominated, 
due to the inefficiency of weak gravitational potentials in 
these systems being able to capture and contain baryonic 
material and then efficiently form stars. Typically these sys- 
tems are strongly disc-dominated, late-type galaxies. 

Intermediate mass systems (11 < log(M v / Mq) < 12) 
generally contain a range of morphological types with sig- 
nificant spheroid and disc components. The stellar content in 
these systems is a strong function of the mass, whereby the 
larger systems become richer in stellar mass fraction reach- 
ing a 'peak' ratio of stellar to dark matter where there is 
a minimum efficiency in the combined effects of supernovae 
and nuclear feedback processes. When we approach the high- 
est mass DM halos (log(M v /M Q ) > 12), feedback through 
nuclear activity generated by an accreting SMBH is able to 
effectively remove the material. This effect becomes stronger 
in the higher mass systems and thus causes the flattening 
and slight downturn of the stellar conversion efficiency. 

Over-plotting model outputs using the 'standard', 'UV 
and 'sophisticated' models, we may easily differentiate be- 
tween outcomes due to the increased sensitivity of the plot- 
ting ranges: We see that the models which include the effects 
of a UV background show significant improvements over the 
entire mass range, and using the two-phase ISM physics and 
adopting the SFR dependent on H2 mass, we are able to re- 
produce the observations to a high accuracy throughout the 
entire mass range due to the relative decrease in SF effi- 
ciency, with the only discrepancies at the lowest masses. It 
is important to note here that due to the sequential build-up 
of matter within halos from high redshifts to low, the im- 
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Figure 10. HI -to-stellar mass ratio as a function of the stellar 
mass. Using derived relationships from §3 we are able to assess 
the relationship between HI and stellar mass, finding that using 
a two-phase ISM and related SFR computations we accurately 
reproduce the observations across the entire observational range, 
unlike the more simplistic frameworks, where we systematically 
convert too much cold gas into stellar material. The shaded region 
represents the 1 — <r observational uncertanty. Model error-bars 
represent poissonian uncertanties due to the finite sample size of 
synthetic galaxies. 

portance of ionizing backgrounds and star formation on low 
mass systems manifests throughout the entire mass range 
at z — (see Somerville et al. 2008 for further discussion 
of this). Under the standard ACDM scenario whereby small 
structures collapse first and form progressively larger sys- 
tems, large structures at z = are therefore significantly in- 
fluenced by small scale processes at higher redshifts, where 
their constituent parts were forming and evolving, thus the 
low mass behavior is of global importance to galaxy forma- 
tion theory. 

4.5 Gas-to-Dark matter properties 

In analogy to the previous section, in Fig. 9 we plot the frac- 
tion of HI to DM halo mass we are able to see that the 
gas fractions are highly sensitive to the different physical 
prescriptions, the scatter is attributed to the fact that the 
cold gas (and thus HI mass) component at any time is con- 
trolled by the competing processes of infall, star formation, 
feedback and recycling. We also find that, due to the uncer- 
tainties and intrinsic dispersion in the relationship between 
the dark matter and gaseous matter, the observational trend 
is not as constrained for the gaseous as for the stellar com- 
ponent. Observationally, the lowest mass DM halos become 
progressively more depleted in HI , due to suppression of in- 
fall and increased efficiency in generating gaseous outflows, 
intermediate mass DM halos (II > (log(M v /M Q ) > 11.5)) 
possess the most HI mass, due to the global maxima in trap- 
ping and containing baryonic material in DM halos (or con- 
versely, a global minima in feedback processes) , and at high 
masses we see a progressive loss of HI abundance in halos, 
signifying an increase in removal of gaseous material through 
cooling inefficiencies and nuclear feedback processes. 

In Fig. 9 we show the outputs for each of the three re- 
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alisations, as can be clearly seen, the only successful repro- 
duction of the majority of the observations is the 'sophisti- 
cated' model, whereby we model both the effects of a two- 
phase ISM and of an ionizing background radiation. There- 
fore modeling the star formation rate depending linearly on 
the surface density of molecular gas in the disc component 
results in a correct conversion of cold gas to stellar mate- 
rial (as can be seen from simultaneous fits to both Fig. 8 
& Fig. 9), using the standard Schmidt-Kennicutt (1998) star 
formation law (and their constrained normalisation) we con- 
vert too much material from gas into stars, resulting in an 
over-production of stellar mass and an underproduction of 
gas mass in all DM halos, whereas the modification outlined 
in Dutton & van den Bosch (2009) allows for the correct 
fraction of virial mass in the HI phase, apart from in the 
lowest mass systems (11 > (log(M v /Mo), where we find ob- 
servation and model predictions in alarming disagreement. 

As with previous works, but to a lesser extent, we find 
that we over-produce the amount of gas in the lowest mass 
DM halos, and more worryingly, we do not reproduce the 
observed down-turn in gaseous mass within the lowest mass 
systems (log(M v /Mq) < 11). This mismatch clearly sug- 
gests that despite our suppression of baryonic material onto 
low mass DM haloes, and improved star formation prescrip- 
tions we are still missing extra factors which dominate at 
M v (z — 0) < 10 n AfQ. We my attribute this effect with our 
lack of description of environmental effects, since the lowest 
mass halos are more probable to be located within larger 
over-densities and therefore subject to external forces (see 
the §5 and Mo et al. 2005 for more details). 

4.6 Stellar-to-Gas properties 

It is suspected that within the lowest mass systems, gas 
collapse, star formation and supernovae feedback result in 
a self-regulated conversion of gas into stars, thus these ef- 
fects can be seen to produce precisely the correct fraction 
of stars to gas in Fig. 10. Within our model the cold gas 
component fluctuates more than any other galaxy compo- 
nent since it is constantly being replenished due to infall and 
stellar recycling, exhausted through star formation and ex- 
pelled through feedback processes. Therefore, directly com- 
paring each galaxies gaseous and stellar properties provides 
a stringent comparison between observation and model. In- 
terestingly we find that plotting these quantities, we repro- 
duce a tight correlation with relatively little scatter. 

We find an overall agreement to the data across the en- 
tire mass range under investigation using the 'sophisticated' 
model, showing that the HI fraction becomes increasingly 
large with decreasing stellar mass, having roughly equal stel- 
lar and HI masses at M s ~ 10 9 Mq increasing to a factor 
of ~ 100 times more HI than stellar mass within systems 
with M B « 1O 7 M0. Conversely, without using the effects of 
a two-phase ISM and UV background radiation, we find a 
general offset in all masses, with an over-efficient conversion 
of cold gas to stars (as may also be seen in Figs. 8 & 9). 

We also find that this result is relatively robust against 
parameter choices; only having a significant dependence on 
the star formation efficiency parameter (which is constrained 
observationally). We may attribute this to the fact that star 
formation is modeled as a function of the surface density of 
H2 in galaxies and within the low mass disc dominated sys- 
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Figure 11. The predicted SFR as function of the stellar mass at 
z = for our three realisations compared to Elbaz et al. 2007. As 
can me seen, the star formation rates in all three models are sim- 
ilar, all showing good agreements in the low/intermediate mass 
systems but models overpredicting the SFR in the highest mass 
systems. The shaded region represents the 1 — a observational un- 
certanty. Model error-bars represent poissonian uncertanties due 
to the finite sample size of synthetic galaxies. 

terns is relatively inefficient (due to the low surface densities 
and relatively high star formation thresholds) resulting in a 
well defined conversion of cold gas to stars and hence, HI to 
stellar mass ratio. This conclusion also has important conse- 
quences when considering the overproduction of both stellar 
and HI masses in the lowest mass DM halos, since we ac- 
curately reproduce the self-regulation properties of galaxies 
(as the correct balance between star formation and feedback 
is required to simultaneously reproduce the correct gas and 
stellar mass budgets in DM haloes). Overproduction in the 
low mass regions within Fig. 8 & Fig. 9 therefore, must be at- 
tributed to the accretion of too much material, a conclusion 
which is confirmed by several other authors (see Mo et al. 
2005 & references therein). 

We would like to plot also the H2 mass counterparts 
for figures 8 & 9, however, due to the unconstrained mass 
functions, as discussed in §4.3, we are not able to construct 
conclusive observational constraints. Secondly, since parti- 
tioning the ISM within the single-phase ISM models is typ- 
ically something relatively ad-hoc, we do not perform this 
analysis here. 

4.7 Star formation properties 

A final useful diagnostic to be used to constrain galaxy 
formation models is the instantaneous star formation rate 
(SFR) which is found to vary significantly within galaxies 
of different stellar mass. In Fig. 11 we compare each model 
realisation with the observational SFR estimates as a func- 
tion of stellar mass by Elbaz et al. 2007, who used a large 
sample of SDSS galaxies with spectroscopic data in order to 
accurately determine the SFR (see Brinchmann et al. 2004). 

We find little difference between the SFR predicted by 
each model realisation as expected since each star formation 
rate prescription is constrained by z = galaxy proper- 
ties. Moreover, we generally find that within the low mass, 
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disc-dominated region [M s < 10 10 Mq) we obtain a good 
agreement between model and observation but we signifi- 
cantly over-predict the star formation rates in the largest 
galaxies. We may attribute this to a lack of AGN quenching 
within the disc component of the largest mass systems at 
late times, which has been applied in several models as a 
rather ad-hoc 'radio-mode' feedback (see Croton et al. 2006, 
Bower et al. 2006). Since, within our models we do not 
have any suppression of the growth of discs around large 
pre-formed bulges, aside from the late transition from the 
spheroid formation epoch to the disc formation epoch. We 
hope to investigate the effects of energetic feedback from 
a formed SMBH-spheroid system on the late properties of 
the disc since we hypothesise that; despite having a negli- 
gible effect on the growth of a pre-formed SMBH, star for- 
mation is prevented within the spheroid component at late 
times since even arbitrarily low accretion rates onto the cen- 
tral SMBH results in energetic feedback capable of heating 
the ISM, however, within our framework, during the quies- 
cent disc growth phase no material is assumed to collapse 
onto the spheroid structure, even with a slight adjustment 
to our model we may allow for some material to collapse 
at late times onto the spheroid-SMBH system, resulting in 
the quenching of star formation in the larger systems and 
naturally generating Seyfert-type galaxies. However, for sim- 
plicity we have neglected this effect within this work, and 
hope to investigate the physical mechanisms capable of gen- 
crating this self-consistently, within a subsequent work. 

4.8 Two component ISM properties 

As a main advance of this model over current SAMs, we 
detail the properties of the ISM by modeling the formation of 
molecular clouds (H2 regions) through pressure arguments 
within the disc component. This enables us to modify the 
star formation law, and thus allows us to gain insights into 
the more detailed gas-properties of normal galaxies under 
the semi-analytical framework. 

Within this work we highlight the importance of mod- 
eling the ISM in two-phases advocating it as a simple, yet 
important advance over current frameworks; whereby the 
formation of H2 regions is determined by the planar pres- 
sure within the gaseous disc structure. This added ingredient 
is important for two main reasons: Firstly the decline in H2 
regions within low mass systems results in a higher fraction 
of cold gas in the form of neutral HI which is thus detectable 
through conventional 21-cm line surveys (see Barnes et al. 
2001), whereas H2 mass estimates prove to be significantly 
more difficult, relying on uncertain conversion factors be- 
tween H2 and CO-lines, therefore, assuming a single con- 
version between total cold gas and HI , as is commonly 
done in SAMs provides inaccurate outputs. Secondly, the 
star formation properties of galaxies have been shown to 
relate explicitly to the detailed properties of the internal 
structure of the ISM (see Krumholz et al. 2009, Gnedin et 
al. 2009, Obreschkow & Rawlings, 2009), and therefore this 
added layer of complexity should now be embedded within 
current SAMs. 

In Fig. 12 we show the fraction of cold gas which is in 
the form of H2 within the disc component as a function 
of total stellar mass (disc and bulge masses). As can be 
seen, there is a tight relationship with little scatter in the 
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Figure 12. The molecular fraction of gas within galaxy disks. 
Showing that in low surface density galaxies, the formation of 
molecular clouds is suppressed, and thus star forming regions 
are diminished, and within the high-mass systems, the spheroid 
component dominates reducing the surface density in disks and 
therefore the molecular gas fraction rapidly declines. We hope to 
compare this result to observational studies as they become avail- 
able. Model error-bars represent poissonian uncertanties due to 
the finite sample size of synthetic galaxies. 



lowest mass systems due to the overall dominance of disc 
components and the molecular fraction decreases steadily 
with decreasing stellar mass, we find a peak in molecular 
fraction corresponding to approximately 80% H2 at M 3 ~ 
1O 11 M0 where the disc mass reaches a maximum. Above this 
mass the spheroid component dominates and the disc surface 
density thus drops, lowering the efficiency of molecular cloud 
formation rapidly. 



5 CONCLUSIONS 

Prompted by several theoretical attempts to model the low- 
mass end of the stellar mass function (Somerville et al. 2008 

6 references therein) and significant observational effort to 
constrain the HI mass function using large surveys of galax- 
ies (Barnes et al. 2001, Zwaan et al. 2005), within this work 
we have developed a physically motivated model in order 
to explain the inefficiencies of low-mass DM halos in trap- 
ping baryonic material and forming stars. Motivated primar- 
ily by the importance of physical descriptions within small 
scale systems which form the building blocks of larger sys- 
tems within a hierarchically clustered Universe, also because 
theoretical models have either neglected, or find significant 
troubles in simultaneously matching both the stellar and HI 
mass functions within the low mass end (see Mo et al. 2005). 

Therefore, within this paper we have utilized well con- 
strained observations of the stellar (Bell et al. 2003a, 2003b) 
and HI (Zwaan et al. 2005) mass functions and employed a 
numerical technique (Shankar et al. 2006) in order to derive 
relationships between galaxy properties and their host DM 
halos. Assuming a one-to-one mapping of these systems, we 
are thus able to make detailed comparisons between models 
and observations. 

In order to interpret these phenomenological relation- 
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ships physically, we develop a cosmologically grounded 
galaxy formation model outlined in C09 but with several 
significant modifications. Under this framework we follow 
the development of baryonic material as it accumulates and 
evolves within growing DM halos and is subject to cooling, 
heating, possible angular momentum losses, star formation 
and recycling, and feedback through both supernovae events 
and the growth of a central supermassive black hole. By 
comparing three model realisations with varying sophisti- 
cation; one using a 'standard' approach, whereby the cos- 
mological baryonic fraction ft = Qb/^im is allowed to ac- 
crete within halos throughout their lifetimes and the star 
formation in discs (and thus the vast majority of low mass 
systems) is given by the Schmidt-Kennicutt (1998) star for- 
mation law, determined by the total cold gas mass, a 'UV 
model, whereby the infall of baryonic material within low 
mass DM halos is suppressed due to the presence of an ion- 
izing UV background (see Eqn.6, 7), and finally a 'sophisti- 
cated' model with a modified infall due to UV radiation and 
a modified star formation law which requires a two-phase 
ISM, comprising of neutral HI and molecular H2 and typi- 
cally results in a lower star formation rate within low mass 
systems at early times. 

These realisations clearly demonstrate that using a 
'standard' approach, the most simple case is not able to si- 
multaneously match both the stellar and HI mass functions, 
significantly over-producing low mass galaxies (sec Figs. 4, 5, 
6, 8 & 9), whereas the use of a 'sophisticated' approach, we 
are able to match observations reasonably accurately to rel- 
atively low masses, finding discrepancies within the lowest 
mass systems. As an additional consequence of the modified 
star formation law in the 'sophisticated' model realisation, 
we naturally partition the cold gaseous ISM into HI and H2 
components, allowing for predictions of both of these quan- 
tities without the need for further assumptions, finding that 
we are able to match both components simultaneously to a 
good accuracy over the entire mass range when comparing 
to the H2 results of (Keres et al. 2003, Obreschkow & Rawl- 
ings, 2009, Fig. 7), and only over-producing the amount of 
HI mass in the lowest mass galaxies (see Zwaan et al. 2005 
& Fig.6). 

Finally, analysing several properties of galaxies against 
their stellar masses, we find that, unlike the simple ap- 
proaches, the HI -to-stellar mass ratio is accurately repro- 
duced using the sophisticated treatment (see Fig. 10), indi- 
cating that the self-regulation of star formation allows for 
the correct conversion of material, is relatively insensitive 
to the infall and supernovae feedback rates within physical 
limits. Comparisons of the star formation rates within these 
models however shows little difference at z = as expected, 
and indicates a secondary problem with our simple physi- 
cal model, over-predicting the SFR in the largest systems 
(with M v > 10 11 M Q ), however, these systems are typically 
spheroid-dominated and are therefore relatively insensitive 
to the details of the disc formation recipes. 

It is clear to assess the limitations and the manifesta- 
tions of this relatively simplistic approach to galaxy forma- 
tion modeling, under our improved framework we still re- 
quire large feedback efficiencies in order to sufficiently sup- 
press star formation within the lowest mass systems: As has 
been studied by Mo et al. 2005, preheating through pre- 
virialised structure formation may further reduce the bary- 



onic infall onto the lowest mass halos, further reducing the 
need for such high SN efficiencies, other environmental ef- 
fects such as tidal stripping and harassment may help to 
further improve the theoretical framework. However, despite 
these further degrees of freedom, we also note that changes 
in the initial mass function towards something more top- 
heavy will naturally have a higher supernovae fraction per 
stellar population, this therefore further reduces the need for 
highly efficient supernovae feedback within these models. In 
order to assess these additional effects we would need to ac- 
count for all the environmental effects associated with galaxy 
evolution, expanding our modeling from a single mass accre- 
tion history to a full merger-tree framework, however, this 
adds a great deal more complication and uncertain physics 
(such as the evolution of satellite structures within DM ha- 
los, the merger rates of galaxies, the outcome of galaxy com- 
ponents in mergers of different ratios), this will be the sub- 
ject of a further analysis in a subsequent paper. 

A second, minor shortcoming of our modeling appears 
to come from the complete separation of spheroid and disc 
growth within the two epochs of DM halo growth, whereby 
the disc and spheroid only share material through disc in- 
stabilities which are generally quite rare at late times due 
to the stabilization generated by the pre-formed spheroid. 
Despite showing that our quasi-monolithic scenario for the 
growth of galaxies to show promising results, we also hope 
to include explicitly the effects of merger events and envi- 
ronmental effects in a future work. 

In conclusion, focusing mainly on the low mass galaxy 
population, adopting several theoretical improvements over 
'standard' SAMs, we are able to simultaneously match both 
the stellar, HI and H2 mass budgets within DM halos, 
and the star formation properties of galaxies within the ob- 
servational ranges. This promising result indicates that at 
present, the 'standard' approach to modeling the low mass 
evolution of galaxies is somewhat over-simplified within cur- 
rent SAMs which only have a single-phase ISM, and, due to 
the hierarchical nature of structure formation, may manifest 
as significant tensions in progressively larger systems within 
a full merger-driven framework (see Somerville et al. 2008 
for a discussion). We therefore advocate the use of more 
sophisticated treatments of the interstellar medium within 
current and future SAMs. 

Assessing the limitations of our framework, we conclude 
that further suppression of infall onto the lowest mass sys- 
tems would allow for a further reduction in the need for 
strong supernovae feedback and should further ease ten- 
sions between models and observations, this could only come 
through environmental effects such as tidal shocks or grav- 
itational pre-heating (Mo et al. 2005), however this effect 
has not been studied in detail through hydrodynamic sim- 
ulations and remains to be fully investigated. By adding a 
channel whereby even small amounts of gaseous material 
may be transferred to the spheroid component during late 
times, small amounts of 'radio mode' AGN activity may be 
triggered, little affecting the spheroid or the SMBH masses, 
but significantly lowering the SFR in the discs, preferentially 
at large masses, naturally resulting in Seyfert-type active 
galaxies and reducing the SFR in these large discs, hopefully 
bringing Fig. 11 into better agreement with observations, we 
also hope to investigate the pan-redshift galaxy population 
under this framework (Cook et al., 2009b submitted). 
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Interestingly however, within this relatively simplistic 
framework we are able to self-consistently reproduce several 
of the key observations, it is therefore clear that mergers, 
to some degree, are not the dominant driver for the global 
evolution of the galaxy population. It will therefore provide 
a useful exercise to mount our physical prescriptions onto a 
full merger-tree DM background which should allow us to 
model environmental effects consistently. The main results 
from this paper indicate however, that using a relatively sim- 
ple framework, we find an reasonable agreement to the stel- 
lar, HI and H2 mass functions of galaxies arising naturally, 
thus, we advocate all current SAMs to begin to incorporate 
two-phase ISM physics into their frameworks. 
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